%##############################################################
%####
%####	Navier-Stokes problem
%####
%##############################################################

% #BEGIN PROBLEM DEFINITION
% #-------------------------------------------------------------
% #Convention used
% # Ny   x  x  x  ... x 
% # .    .
% # .    .
% # .    . 
% # 3    x  x  x  ... x  
% # 2    x  x  x  ... x  
% # 1    x  x  x  ... x  
% #
% # j/i  1  2  3  ... Nx

tic

%% Convergence parameters--------------------------------------------------
convTol = 10^-4;
%--------------------------------------------------------------------------
L = 1.0;
global N;
N = 11;

visc = 0.01;
h = 1./N

rho = 1.0;

dimp = 20*N^2;
dimu = dimp - 2*N;
dimv = dimp - 12*N;
%# initial data
u = zeros(dimu,1);
v = zeros(dimv,1);
p = zeros(dimp, 1);
%# 
unovo = zeros(dimu,1);
vnovo = zeros(dimv,1);
%#


%#temporal integration
dt = 0.001;
T0 = 0.;
T1 = 80.;

if( min(h^2/(4*visc), 2*visc) < dt) 
	fprintf ('\n\t\t ^^^*  No convergencia  ^^^*\n\n')
end


NumTime = (T1 - T0)/dt;

%%% Salvar solução cada 0.1 segundo
plot_t = 0.1/dt;
usol = zeros(dimu, T1/0.1 + 1);
vsol = zeros(dimv, T1/0.1 + 1);
psol = zeros(dimp, T1/0.1 + 1);

usol(:,1) = u;
vsol(:,1) = v;
psol(:,1) = p;
temp = 2;


%## Matrix for pressure
Ap = zeros(dimp, dimp);
bp = zeros(dimp, 1);

for i = 1: 12*N
	if(i < 4*N + 1) jmin = N + 1;
	else jmin = 1;
    end
	for j = jmin: 2*N
		centro = 4;
		if( j > jmin) Ap(pij(i,j),pij(i,j-1)) = 1.0; else centro = centro - 1; end
		if( j < 2*N ) Ap(pij(i,j),pij(i,j+1)) = 1.0; else centro = centro - 1; end
		if( i > 1)    Ap(pij(i,j),pij(i-1,j)) = 1.0; else centro = centro - 1; end
		if( i < 12*N) Ap(pij(i,j),pij(i+1,j)) = 1.0; else centro = centro + 1; end
		Ap(pij(i,j),pij(i,j)) = -centro;
	end
end
Ap = Ap/(h^2);
Ap = sparse(Ap);

%#### Temporal loop
for tt = 1: NumTime
	fprintf('     SOLVING.. %d iteration\n', tt);
	
	%####   Predictor step   #######
	for i = 1: 12*N
		if(i < 4*N + 1) jmin= N + 1;
		else jmin = 1;
	    end
		for j = jmin: 2*N
			if(i < 12*N) 
				upto = u(uij(i,j),1); 
				yy = 1 + 0.5*h + h*(j - jmin);
				
				if(i < 12*N - 1)	uri = u(uij(i+1,j),1); else uri = upto; end
				if(j == 2*N) uup = -upto; else uup = u(uij(i,j+1),1); end
					
				if(i == 1) ule = 4*(yy - 1)*(2 - yy);
				elseif(i == 4*N + 1 && j < N + 1) ule = 0.0;
				else ule = u(uij(i-1,j),1); 
				end
				
				if(j == jmin) udo = -upto; 
				else udo = u(uij(i,j-1),1);
				end
						
				if (i == 1) udi = 4*(yy+h - 1)*(2 - yy-h);
				elseif ( i == 4*N + 1 && j < N ) udi = 0.0;
				elseif ( j == 2*N) udi = -u(uij(i-1,j),1);
				else udi = u(uij(i-1,j+1),1);
				end
			end
			if(j < 2*N) 
				vpto = v(vij(i,j),1); 
				
				if( i < 12*N)		vri = v(vij(i+1,j),1); else vri = -vpto; end
				if( j == 2*N - 1)   vup = 0.0; else vup = v(vij(i,j+1),1); end
				
				if( i == 1) vle = -vpto;
				elseif(i == 4*N + 1 && j < N + 1) vle = -vpto;
				else vle = v(vij(i-1,j),1);
				end
				
				if( j == jmin ) vdo = 0.0; else vdo = v(vij(i,j-1),1); end
				
				if(j == jmin) vdi = 0.0;
				elseif(i == 12*N) vdi = -v(vij(i, j-1),1);
				else vdi = v(vij(i+1,j-1),1);
				end
			end
			
			
			Ax = ((uri + upto)^2 - (upto + ule)^2 + (uup + upto)*(vri + vpto) - (upto + udo)*(vdi + vdo))/(4*h);
			Ay = ((vup + vpto)^2 - (vpto + vdo)^2 + (uup + upto)*(vri + vpto) - (udi + ule)*(vpto + vle))/(4*h);
			
			Dx = (uri + ule + uup + udo - 4*upto)/(h^2);
			Dy = (vri + vle + vup + vdo - 4*vpto)/(h^2);
			
			if(i < 12*N)
				unovo(uij(i,j),1) = upto + dt*(-Ax + visc*Dx);
			end
			if(j < 2*N)
				vnovo(vij(i,j),1) = vpto + dt*(-Ay + visc*Dy);
			end
			
		end
	end
	
	%####   Pressure Poisson equation  #######
	for i = 1: 12*N
		if(i < 4*N + 1) jmin = N + 1;
		else jmin = 1;
	    end
		for j = jmin: 2*N
			if(i == 12*N) uri = unovo(uij(i-1,j),1); else uri = unovo(uij(i,j),1); end
			
			if(i == 1) yy = 1 + 0.5*h + h*(j - jmin); ule = 4*(yy-1)*(2 - yy);
			elseif(i == 4*N + 1 && j < N + 1) ule = 0.0;
			else ule = unovo(uij(i-1,j)); end
			
			if(j == 2*N) vup = 0.0; else vup = vnovo(vij(i,j)); end
			
			if(j == jmin) vdo = 0.0; else vdo = vnovo(vij(i,j-1)); end
			
			bp(pij(i,j),1) = rho*(uri - ule + vup - vdo)/(dt*h);
		end
	end
	%#p = gmres(Ap, bp, rtol = 1e-2);
	p = Ap\bp;
	
	%####   Predictor step  #######
	for i = 1: 12*N
		if(i < 4*N + 1) jmin = N + 1;
		else jmin = 1;
	    end
		for j = jmin: 2*N
			ppto = p(pij(i,j));
			if(i < 12*N)
				upto = unovo(uij(i,j),1);
				pri = p(pij(i+1,j),1);
				u(uij(i,j),1) = upto - dt*(pri - ppto)/(rho*h);
			end
			if(j < 2*N)
				vpto = vnovo(vij(i,j),1);
				pup = p(pij(i,j+1),1);
				v(vij(i,j),1) = vpto - dt*(pup - ppto)/(rho*h);
			end
		end
    end
    
    if(mod(tt,plot_t) == 0) 
        usol(:,temp) = u(:,1);
        vsol(:,temp) = v(:,1);
        psol(:,temp) = p(:,1);
        temp = temp + 1;
    end
    
    %% convergence check---------------------------------------------------
    if (tt > 1)
        uDiffSum = 0;
        for i = 1 : dimu;
            uDiff = abs(u(i) - um1(i));
            if (u(i) ~= 0)
                uDiff = abs( (u(i) - um1(i)) / u(i) );
            end;
            uDiffSum =  uDiffSum + uDiff;
        end;        
        vDiffSum = 0;
        for i = 1 : dimv;
            vDiff = abs(v(i) - vm1(i));
            if (v(i) ~= 0)
                vDiff = abs( (v(i) - vm1(i)) ./ v(i) );
            end;
            vDiffSum =  vDiffSum + vDiff;
        end;        
        err = uDiffSum/dimu + vDiffSum/dimv;
        fprintf('Convergence error: %f \n', err);
        if (err <= convTol)
            fprintf('Stationary state reached at iteration %d, convergence error: %f', tt, err);
            convFlag = true;
            break;
        end;
    end;
    % the previous time step velocities.
    um1 = u;
    vm1 = v;
    %----------------------------------------------------------------------
end

toc

%##############   BEGIN POST-PROCESSING STEP

%### dom 1
x1 = zeros(4*N + 1, 1);
y1 = zeros(N + 1, 1);
x1(1) = 0.0;
y1(1) = 1.0;
for( i = 2: 4*N + 1 ) x1(i) = x1(i - 1) + h; end
for( i = 2: N + 1)    y1(i) = y1(i - 1) + h; end

u_dom1 = zeros((4*N + 1)*(N + 1),1);
v_dom1 = zeros((4*N + 1)*(N + 1),1);
p_dom1 = zeros((4*N + 1)*(N + 1),1);

x2 = zeros(8*N + 1, 1);
y2 = zeros(2*N + 1, 1);
x2(1) = 4.0;
y2(1) = 0.0;
for( i = 2: 8*N + 1 ) x2(i) = x2(i - 1) + h; end
for( i = 2: 2*N + 1)  y2(i) = y2(i - 1) + h; end

u_dom2 = zeros((8*N + 1)*(2*N + 1),1);
v_dom2 = zeros((8*N + 1)*(2*N + 1),1);
p_dom2 = zeros((8*N + 1)*(2*N + 1),1);

[X1, Y1] = meshgrid(x1, y1);
[X2, Y2] = meshgrid(x2, y2);
figure()

for(tpo = 1: T1/0.1 + 1)
    u(:,1) = usol(:, tpo);
    v(:,1) = vsol(:, tpo);
    p(:,1) = psol(:, tpo);
    kk = 0;
    for( i = 1: 4*N + 1)
        for (j = N + 1: 2*N + 1)
            kk = kk + 1;
            if(j == N + 1)
                u_dom1(kk) = 0.0;
                v_dom1(kk) = 0.0;
                if(i == 1) p_dom1(kk) = p(pij(i,j));
                else p_dom1(kk) = 0.5*(p(pij(i-1,j)) + p(pij(i,j)));
                end
            elseif(j == 2*N + 1)
                u_dom1(kk) = 0.0;
                v_dom1(kk) = 0.0;
                if(i == 1) p_dom1(kk) = p(pij(i,2*N));
                else p_dom1(kk) = 0.5*(p(pij(i-1,j-1)) + p(pij(i,j-1)));
                end
            elseif(i == 1)
                yy = 1 + h*(j - N - 1);
                u_dom1(kk) = 4*(yy-1)*(2 - yy);
                v_dom1(kk) = 0.0;
                p_dom1(kk) = 0.5*(p(pij(i,j-1)) + p(pij(i, j)));
            else
                u_dom1(kk) = 0.5*(u(uij(i-1,j-1)) + u(uij(i-1,j)));
                v_dom1(kk) = 0.5*(v(vij(i-1,j-1)) + v(vij(i,j-1)));
                p_dom1(kk) = 0.25*(p(pij(i-1,j-1)) + p(pij(i,j-1)) + p(pij(i-1,j)) + p(pij(i,j)));
            end
        end
    end

    kk = 0;
    for( i = 4*N + 1: 12*N + 1)
        for (j = 1: 2*N + 1)
            kk = kk + 1;
            if(j == 1) 
                u_dom2(kk) = 0.0;
                v_dom2(kk) = 0.0;
                if(i == 4*N + 1) p_dom2(kk) = p(pij(i,j));
                elseif(i == 12*N + 1) p_dom2(kk) = 0;
                else p_dom2(kk) = 0.5*(p(pij(i-1,j)) + p(pij(i,j)));
                end
            elseif(j == 2*N + 1)
                u_dom2(kk) = 0.0;
                v_dom2(kk) = 0.0;
                if(i == 1) p_dom2(kk) = p(pij(i,2*N));
                elseif(i == 12*N + 1) p_dom2(kk) = 0;
                else p_dom2(kk) = 0.5*(p(pij(i-1,j-1)) + p(pij(i,j-1)));
                end
            elseif(i == 4*N + 1 && j < N + 2)
                u_dom2(kk) = 0.0;
                v_dom2(kk) = 0.0;
                p_dom2(kk) = 0.5*(p(pij(i,j)) + p(pij(i,j-1)));
            elseif(i == 12*N + 1)
                v_dom2(kk) = 0.0;
                u_dom2(kk) = u(uij(i-2,j-1));
                p_dom2(kk) = 0;
            else
                u_dom2(kk) = 0.5*(u(uij(i-1,j-1)) + u(uij(i-1,j)));
                v_dom2(kk) = 0.5*(v(vij(i-1,j-1)) + v(vij(i,j-1)));
                p_dom2(kk) = 0.25*(p(pij(i-1,j-1)) + p(pij(i,j-1)) + p(pij(i-1,j)) + p(pij(i,j)));
            end
        end
    end


    uu1 = reshape(u_dom1, N + 1, 4*N + 1);
    vv1 = reshape(v_dom1, N + 1, 4*N + 1);
    pp1 = reshape(p_dom1, N + 1, 4*N + 1);

    uu2 = reshape(u_dom2, 2*N + 1, 8*N + 1);
    vv2 = reshape(v_dom2, 2*N + 1, 8*N + 1);
    pp2 = reshape(p_dom2, 2*N + 1, 8*N + 1);

    tit = ['Tempo = ', num2str((tpo-1)*0.1)]; 
    
    if(tpo < 10) tit2 = ['NS_000',num2str(tpo),'.png'];
    elseif(tpo < 1000) tit2 = ['NS_00',num2str(tpo),'.png'];
    else tit2 = ['NS_0',num2str(tpo),'.png'];
    end

    subplot(3,1,1)
    surf(X1, Y1, uu1)
    %pcolor(X1, Y1, uu1)
    hold on
    surf(X2, Y2, uu2)
    % pcolor(X2, Y2, uu2)
    hold off
    shading interp
    view(0, 90)
    ylabel('Horizontal velocity')
    if(tpo == T1/0.1 + 1) colorbar('location', 'North'); end
    title(tit)

    subplot(3,1,2)
    surf(X1, Y1, vv1)
    hold on
    surf(X2, Y2, vv2)
    hold off
    shading interp
    view(0, 90)
    ylabel('Vertical velocity')
    if(tpo == T1/0.1 + 1) colorbar('location', 'North'); end
    %colorbar();
    
    subplot(3,1,3)
    surf(X1, Y1, pp1)
    hold on
    surf(X2, Y2, pp2)
    hold off
    shading interp
    view(0, 90)
    ylabel('Pressure')
    if(tpo == T1/0.1 + 1) colorbar('location', 'North'); end
    %colorbar();
    
    
    saveas(gcf, tit2);
    pause(0.02)
end